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ABSTRACT 

A numerical code for solving various Lya radiative transfer (RT) problems is presented. The code is 
suitable for an arbitrary, three-dimensional distribution of Lya emissivity, gas temperature, density, 
and velocity field. Capable of handling Lya RT in an adaptively refined grid-based structure, it enables 
detailed investigation of the effects of dumpiness of the interstellar (or intcrgalactic) medium. The 
code is tested against various geometrically and physically idealized configurations for which analytical 
solutions exist, and subsequently applied to three different simulated high-resolution "Lyman- break 
galaxies", extracted from high-resolution cosmological simulations at redshift z = 3.6. Proper treat- 
ment of the Lya scattering reveals a diversity of surface brightness (SB) and line profiles. Specifically, 
for a given galaxy the maximum observed SB can vary by an order of magnitude, and the total flux 
by a factor of 3-6, depending on the viewing angle. This may provide an explanation for differences 
in observed properties of high- redshift galaxies, and in particular a possible physical link between 
Lyman-break galaxies and regular Lya emitters. 

Subject headings: galaxies: formation — galaxies: evolution — galaxies: fundamental parameters 
(classification) — radiative transfer — scattering — line: formation — line: profiles 



1. INTRODUCTION 

The significance of the Lya emission line as a probe 
of the high-redshift Un iverse has long been establ ished. 
In their classic paper, iPartridge fc Peebles] (|l967t ) sug- 
gested how detection of young galaxies would be fea- 
sible using the Lya line. Nevertheless, for almost 
three decades only a few Lya emitters (LAE s ) wer e 
discovered (see, e.g., iDiorgovski fc Thompson! Il992f ). 
Several theories were proposed to explain the many 
null res ults, e.g. suppression of the Lya line due to 
metals (iMeier fc Terlevichl 1 1 9 8 If) . absorption by dust 
(jHartmann et al.lll988f) and lower-than-ex pected forma- 
tion of massive stars TValls GabaudlfT993]) . 

However, as surveys eventually were able to go deeper, 
and as searching wide regions on the sky became feasi- 
ble, large numbers of high-redshift, star-forming galaxies 
were discovered. Notable s urveys include the La rge Area 
Lyman Alpha survey (e.g.. iRhoads et aT]|2000D and th e 
Subaru Deep Field survey (e.g.. iTanieuchi et alJ l2005h . 
Currently, rcdshifts of LAEs up to z ~ 7 (|Ive et al 120071) 
have been reached, while galaxies detec ted via z-band 
dropo ut observations have reached z ~ 8 (jBouwens et alJ 
2004), but in time deeper observations will be real- 
ized, with the advent of, e.g., the Ultra- VISTA project 
(to be launched u ltimo this year, reaching z = 8.8; 
iDunlop et"aI1 120071 ) and the Jame s Webb Space Tele- 
scope (to be launched in 2013; e.g.. lGardnerll2006l ). 

Besides contributing to our understanding of the over- 
all structure and evolution of the Universe, much in- 
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sight into the properties and formation of the galaxies 
themselves, the fundamental building blocks of our Uni- 
verse, has now been gained from these high-redshift sur- 
veys. Thus, num erous characteri stics of the LAEs, such 
as their density (iHu et al.1ll998[ ). clustering properties 
(lOuchi et al.l 120031 ). and luminosity function (|Hu et alJ 
2003), have been subject to investigation. 

Radiative transfer (RT) is playing an increasingly im- 
portant role in numerical astrophysics and cosmology. 
This is particularly true in the case of Lya. Due to 
the resonant scattering nature of Lya radiation, and due 
to the fact that neutral hydrogen is abundantly present 
in the interstellar medium (ISM) and the intergalactic 
medium (IGM), the history of a Lya photon, from the 
moment it is created in a high-redshift galaxy to the 
time it is captured by a telescope, is not trivial. In gen- 
eral, analytical solutions of RT problems are obtainable 
only in very idealized cases. By far, most of the work 
done on the subject has been concerned with the emerg- 
ing spectrum from an isothermal, homogeneous m edium 
of p l ane-parallel or spher i cal symmetry (e.g. lAuerl 
1965t lAverv fc House! [ l968t iPanangia fc Ranieril 11973 ' 



Ahn et all 120011 120021 IZheng fc Miralda-Escuddl2002¥ 
Some allow for isotropi c velo c ities (e.g. ICaroff et all 
1971 iNatta fc Beckwithl Il986l lLoeb fc Rvbickil 119991 : 



Diikstra et all 120061). and some include simple mod 
els for dust (e.g. iBonilha et all fT97j; |Ahja_gtjai1 l2000t 
lHansen fc Obll2006l IVerhamme e t al. 2006, 2Q0l). How- 



ever, even though the results of this work have improved 
tremendously our knowledge of many physical processes, 
they do not capture the complexity and diversity of re- 
alistic, astrophysical situations where velocities can be 
quite chaotic, and densities and temperatures can vary 
by many orders of magnitude over relatively small dis- 
tances. 

To this end, a few codes with varying aims have been 
constructed, capable of performing realistic RT for ar- 
bitrary distributions of source Lya emission, neutral 
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hydrogen density, temperature and velocity resulting 
from cosmological simulations so as to yield the spec- 
trum and the spatial distribution of the escaping pho- 
tons (|Cantalupo et alj|2 005: Tasitsiomi 2006a; Kollmeier 
|2006| ). Also IVerhamme et al.l ()2006l ) have presented a 
similar code, although it has not been applied to cosmo- 
logical simulations. 

Although the work carried out in this paper is largely 
inspired by previous studies, it improves on earlier works 
in several differ ent ways: most imp ortantly, as in the case 
of the code of iTasitsiomil (|2006al ). our code is capable 
of working with physical data on an adaptively refined 
mesh, as opposed to a regular grid. Since gas clumping 
affects the photon escape probability, very high resolu- 
tion is desired. Our study will be restricted to quite 
evolved galaxies, on the kpc scale. Furthermore, in ad- 
dition to studying the emergent spectrum and surface 
brightness distribution, we will investigate the effect of 
viewing the systems from different angles. 

The rest of the paper is organized as follows: in fJ2J a 
basic theory of Lya RT is briefly described. The equa- 
tions presented in this section will serve as a basis for 
understanding and testing the code. A detailed descrip- 
tion of the principles of the RT code is given in SJ31 and in 
2] tests of the code against various analytical solutions 
are presented. A semianalytical acceleration scheme is 
derived in SJ3J and in S^thc code is applied to three differ- 
ent galaxies extracted from high-resolution cosmological 
simulations. Finally, a discussion of the results is given 
in d 

2. RESONANT SCATTERING RADIATIVE TRANSFER 

The first attempts to predict the diffusion of Lya were 
made under the as sumption of coherent scattering in 
the o bservers frame (jAmbarzumi an 1932; Cha ndrasekharl 
Il935t) . The probability of interaction between a photon 
and an atom at rest with respect to the reference frame 
in which the frequency v of the photon is measured is 
described by the line profile (f>(v)\ this was known to be 
given by the sharply peaked natural (Lorentzian) line 
profile C of width A^l = 9.936 x 10 7 Hz around the line 
center frequency vq = 2.466 x 10 15 Hz. 

Several physical quantities that may or may not be 
directly observable have been the subject of interest, 
e.g. the average number (A scat ) of scatterings required 
to escape the medium (to determine the probability of a 
photon being destroyed by dust grains, or by collisions of 
the scattering atom with other atoms while being in the 
excited state) and the shape of the emergent spectrum. 
Due to the complexity of the problem, the physical 
configurations investigated have traditionally been 
restricted to problems in which photons are emitted in 
the center of a homogeneous, isothermal cloud which is 
either spherically symmetric, or infinite in two directions 
and has a finite extension in the third direction (a 
plane-parallel "slab" ) . Denoting by To the optical depth 
for a photon in the line center from the initial point of 
emission to the edge of the gaseous cloud, from pure 
random walk considerations one would naively infer 
(A^cat) ~ Tq. Accordingly, the medium would not have 
to be very opaque in order for the destruction processes 
of Lya to become significant. 



2.1. Thermal Broadening of the Line 

iHenvevI (|1940D and ISpitzerl (|1944f ) acknowledged the 
fact that scattered photons undergo a change in fre- 
quency due to thermal Doppler broadening of the scat- 
tering atoms. In the following, to simplify notation the 
frequency of the photon is parametrized through x = 
(y — i>o)/Avd, where Aisd — (fth/c)^ is the width of the 
Doppler (Gaussian) profile, with w t h = (SfcsT/mH) 1 ^ 2 
being the thermal atom velocity dispersion (times \f2) 
and the rest of the variables having their usual meaning. 
In terms of these quantities, with <\>{y)dv = <p(x)dx the 
normalized thermal line profile is 

Q{x) = ^=e-*\ (1) 

while the natural line profile is 

C{x) = -^ r5 , (2) 
TT x z + a z 

where a = Avl /2Avd is the relative line width. The 
resulting (Voigt) profile is a convolution between the two, 
but due to the smallncss of a, the center of the profile is 
entirely dominated by Q. 

Relying on these considerations. IZanstral (j 19491 . fl95ll ) 
argued that, in each scattering, the frequency of the Lya 
photon would undergo complete redistribution over the 
Doppler line profile, i.e. there is no correlation between 
the frequency Xi of the incoming and xj of the outgo- 
ing photon, and the probability that x < Xf < x + dx 
is 4>{xf)dx. In this picture, the photon still executes 
a random walk, but at each scattering there is a small 
possibility that it will be redistributed so far into the 
wing as to render the medium optically thin and thus 
allow escape. This reduces (A scat ) significantly, and the 
result was later verified n umerically for inte rmediate op- 
tical depths (r - 10 4 ) bv lKoelbloedl (fl95l . 

Stil l based on the assumption of isotropic scattering, 
lUnnol (|1952af bl) calculated an "exact redistribution" for- 
mula q(xi,Xf), giving the probability distribution of x f 
as a function of Xi. With this result. lOsterbrockl (|1962f ) 
found that in the wings, the rms frequency shift (Aa;) rms 
per scattering is 

(Ax) rma = 1, (3) 

and the mean shift (Ax) per scattering is 

(Ax) = -l/M, (4) 

i.e. there is a tendency to drift toward the line cen- 
ter. Thus, a photon at frequency x S> 1 will execute 
a nearly random walk in frequency, returning to the core 
in A'scat.ret. ~ x 2 scatterings. 

From g(xi, z/), IO~sterbrockl (|1962f ) found that (A scat ) cx 
To for moderate optical depths. However, he argued that 
for some limiting large optical depth — which he was 
not able to calculate due to the lack of "sufficiently large 
digital computers" — Xf is so large that the photon will 
execute a random walk also in real space, whence in this 
case (A/scat) oc t$. 

None t heless , applying the method of lFeautrierl (|1964f) . 
I Adams! (|1972f) found numerically that also for extremely 
large optical depths (t up to 10 8 ), (A scat ) oc t . Al- 
though he could not prove it rigorously, he was able to 
give a heuristic argument on physical grounds for this 
behavior. 
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2.2. Neufeld Solution 

The resu l t wa s proven the subs e quent year by 
iHarringtonl (|1973f ): inspired by lUnnd (|1955[ ). utilizing 
the Eddington approximation — which implies that the 
radiation field is everywhere nearly isotropic, but with a 
small net outward flow — and expanding the rcdistribu- 
tion fu nction as formulat ed bv lHummerl (|l962j) to second 
order, IHarringtonl (|1973l ) obtained a diffusion equation 
for the angular averaged intensity J(r, x) within a (non- 
absorbing) slab of extremely large optical depths (de- 
fined 4 as or > 10 3 /V^, or r > 1.2 x 10 6 for T = 10 4 
K)- 

With the photons emitted isotropically from a central 
source emitting 1 photon per unit time, i.e. 1/4tt photons 
per unit time per steradian, an initial frequency x nl] = 0, 
and scatterings assumed to be do minated by isotropic 
wing scatterings, IHarringtonl (|1973[) obtai n ed an expres- 
sion for the emergent spectrum. " INeufeldl (fl99f)h gave a 
more general solution to the problem, allowing for the 
destruction of photons and the injection at any initial 
optical depth r in the slab, with arbitrary initial fre- 
quency. For centrally 5 (at r = 0) emitted radiation in a 
nonabsorbing medium, the solution at the surface, i.e. at 
r = ±t , is 

j(± T x ) = y^_^ 1 

T °' X> 24 ^ar cosh [yfi&Jte (x 3 - xf nj )/ar ] ' 

With perhaps some injustice, we will refer to Eq. ([5]) as 
the "Neufeld solution", even when Xj n \ = 0, in which 
case it reduces to the result of IHarringtonl ([1973). The 
profile is normalized to 1/4tt and exhibits two bumps, 
symmetrically centered on x = and drifting further 
apart for increasing otq. Note that it solely depends 
on the product aro, and that the physical size of the 
gaseous system does not enter the equation. A higher 
density is compensated for by a higher temperature, since 
aro oc (A^ 1 )(?ihiAi/^ 1 ) oc n^/T (at a given size). 
The physical explanation for this is that the denser the 
medium is, the further into the wing the photons have to 
drift. Meanwhile, a higher temperature — and a result- 
ing higher velocity dispersion of the atoms — will make 
the medium less opaque to radiation, since this means 
fewer atoms with a velocity matching the frequency of 
the photons. 

Setting dJ/dx = , IHarringtonl (|1973l ) showed that 
the emergent spectrum has its maximum at 

x m = ±1.066(aT ) 1 / 3 , (6) 

while the average number of scatterings that a photon 
undergoes before escaping the slab was shown to be 

(A scat ) = 1.6127b. (7) 

Except for nume rical factors of order unity, 
iDiikstra et all (|2006l ) derived similar expressions 

4 Note that in [Harringtonj's papers, as well as most coeval au- 
thors', the optical depth at frequency x is defined as t x = tq<P(x), 
whereas in our definition t x = roH(a, x). Since H (a, x) = y/n<j>(x), 
this implies that THarrington = %/7TTus. This definition has been 
chosen to follow more recent studies. 

5 INeufeldl assumed that the photons are emitted from a thin 
layer inside the slab, parallel to the surface. However, for reasons 
of symmetry, we may as well assume that they are emitted from a 
single point. 



for the emergent spectrum and the number of scat- 
terings for photons escaping a static, isothermal, 
homogeneous sphere of gas. Furthermore, the spectrum 
for an isotropically expanding (as in Hubble flow) or 
contracting (as in a gravitational collapse) medium, but 
with no therma l motio n, was examined analytically by 
lLoeb fc Rvbickil |l£99). Evidently, all of the configura- 
tions considered so far are highly idealized compared to 
realistic, astrophysical situations and for more general 
geometries and velocities, analytic solutions arc not 
obtainable. Nevertheless, they provide valuable and at 
least qualitative insight into the characteristics of young 
galaxies, the ISM and IGM, Hi envelopes surrounding 
hot stars, etc. Moreover, they offer direct means of 
testing numerical methods (see 0. 

3. THE CODE 

The transfer of the Lya photons is conducted us- 
ing the 3D adaptive mesh refinement (AMR) Monte 
Carlo code MoCaLaTA. Except for a few numerical 
improvements, in particular the acceleration scheme de- 
scribed in j}5l the code res e mbles the one presented in 
lLaursen fc Sommer-Larsenl ()2007l ). with one significant 
advance: it is now capable of assuming an adaptively re- 
fined mesh, to an arbitrary level of refinement. This al- 
lows for the opportunity to study the effect of the dumpi- 
ness of the ISM on the radiative transfer in great detail. 

The principles of the code wer e briefly explained in 
lLaursen fc Sommer-Larsenl ()2007l ). In the following we 
give a more elaborate description of how the RT is real- 
ized. 

3.1. Lya Emission 

The physical volume of interest is discretized on a base 
grid, typically consisting of 128 3 cells. Cells may be sub- 
divided into eight subcclls which, in turn, may be further 
refined. The refinement criterion is usually taken to be 
density, but can in principle be any condition, e.g. den- 
sity gradient, velocity, etc. If the underlying cosmological 
simulation is particle based, as is the case in the present 
study, the physical parameters of interest are first inter- 
polated onto the grid. 

Each cell contains information about the Lya luminos- 
ity LLy Q and the density nm of neutral hydrogen, as well 
as the temperature T and the three-dimensional peculiar 
velocity field Vbuik of the gas elements. The ratio of L\^ ya 
of a given cell to the total luminosity L to t of all cells de- 
termines the probability of a photon being emitted from 
that particular cell. The initial position Xj of the pho- 
ton is a random location in the cell. In the reference 
frame of the emitting atom, the photon is injected with 
a frequency x na ±, given by the distribution C(x) (Eq. [2]). 
The atom, in turn, has a velocity v at0 m in the reference 
frame of the gas element drawn from a thermal profile of 
Dopplcr width A^d- Measuring atom velocities in terms 
of Dopplcr widths, u = v atom /tith, each component m is 
then distributed according to Q(ui), given by Eq. ([I]). 

The initial direction of the photon follows an 
isotropic probability distribution. To first order in v/c, 
this is true in all relevant reference frames, and a Lorentz 
transformation to the reference frame of the gas element 
then yields the initial frequency xi = x nat + u • fij . 

For photons emitted in the dense, star-forming regions, 
it makes no difference whether x\ is calculated in the 
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above manner or simply set equal to zero. However, when 
studying large volumes of space, a nonvanishing fraction 
of the Lya photons may be produced through cooling 
radiation, which also takes place well away from the star- 
forming regions of the galaxy. In these environments, 
whereas the probability of a photon with x = escaping 
is still extremely small, being injected one or two Dopplcr 
widths away from line center may allow the photon to 
escape. 

3.2. Propagation of the Radiation 

The optical depth r covered by the photon before it is 
scattered is governed by the probability density function 
-P(t) = e~ r , and after initial emission and all subsequent 
scatterings, a random value of r is drawn from P(t). 
This optical depth is converted into a physical distance 
r = t J Tt\\\O x . In the reference frame of the gas, the cross- 
section o~ x of the atom responsible for the scattering event 
is given by the Voigt profile, resulting in 

^ = /i2 H(a,x), (8) 

where /12 = 0.4162 is the Lya oscillator strength, m e is 
the mass of the electron, and 

/+OO — y 2 

is the Voigt function. This function can be approxi- 
mated by a Gaussian in the core and a power law in the 
wing. However, in the transition domain between core 
and wing, either approxi mation is poor. H ence, instead 
we use the analytical fit (|Tasitsi omi 20Q6a|) 

H(a,x) =q^ + e~ x \ (10) 

where 

f for C < 

q = 1 (1 + i) ^TTyn(C) for C >0, (") 

with C = (x 2 - 0.855)/(x 2 + 3.42) and 11(0 = 5.674C 4 - 
9.207C 3 + 4.421C 2 + 0.11170 This is an excellent approx- 
imation for all frequencies at temperatures above 2 K. 

All physical parameters entering the equations above 
are of course given by the cell in which the photon is 
presently located; the host cell. The new position is then 
xj = Xj + rfij. However, since in general the physical 
conditions vary from cell to cell, if xy is outside the host 
cell, the photon is placed at the point x cut of intersection 
with the face of the cell and the above calculation is 
redone with the parameters of the new cell. Part of the 
originally assigned value of r has already been "spent" , 
so the remaining optical depth to be traveled is now 

T = Torig. - I Xcut - X i|( n Hr0-s) prevxeir (12) 

The frequency of the photon is Lorentz transformed to 
the bulk velocity of the new cell. 

In contrast to a regular grid, in an AMR grid a given 
cell will not in general have a unique neighbor. The cells 
are structured in a nested grid, where a refined cell is the 
"parent" of eight "child" cells which, in turn, may or may 
not be refined. The new host cell of the photon is then 
determined by walking up and down the hierarchical tree 
structure. 




Fig. 1. — Illustration of the mechanism responsible for the fre- 
quency shift of a scattered photon. In the reference frame of an 
external observer (left), a photon blueward of the line center (blue 
solid) is scattered by an atom receding in such manner that the 
component un of its velocity u along the direction of the photon 
matches closely the frequency x. In the reference frame of the atom 
(right), the photon then seems close to the line center (green). Ex- 
cept for a minute recoil effect, the photon leaves the atom with the 
same frequency. However, to the external observer, if the photon 
is scattered in a direction opposite the atom's motion (red), it will 
be redshifted. Only if by chance it is scattered in exactly the orig- 
inal direction (dotted blue), its frequency remains unaltered. For 
| a: | S> 1, the number of atoms with «|| ~ x is so small that the 
photon is most likely to be scattered by a low-velocity atom. In 
this case, no matter in which direction the photon is scattered the 
motion of the atom will not shift the frequency significantly. 

3.3. Scattering 

When the initially assigned r is "used up" , the photon 
is scattered. It must be emphasized that the discussed 
broadening of the line and the corresponding uncertainty 
in energy does not imply that a photon of a given en- 
ergy can be absorbed, and subsequently re-emitted with 
a different energy. Indeed, this would be possible had 
the energy of the ground state been associated with an 
uncertainty in energy as well. However, since the life- 
time of this state is effectively infinite, its energy is well- 
defined. Except for a small recoil effect, the scattering 
is coherent in the reference frame of the atom. How- 
ever, to an external observer the nonzero velocity of the 
scattering atom will, in general, add a frequency shift to 
the photon. Figure [T] shows a qualitative interpretation 
of how the Doppler shift arises. Since the frequency de- 
termines the opacity of the gas, the exact value of the 
velocity is important. In the directions perpendicular to 
fi;, the velocities uj_i,2 will follow a Gaussian distribu- 
tion G(u±i2)- However, due to the resonance nature of 
the scattering event, the velocity it|| parallel to ft; de- 
pends on x. Thus, the probability distribution £/(un) 
must be convolved with the probability C(x — uu) of the 
atom being able to scatter the photon. For small val- 
ues of being scattered by an atom with itn = x is 
highly favored. For large values of |x| the abundance 

of these atoms decrease as e~ x , so that scattering by 
"slow" atoms becomes more probable, even though in 
reference frame of these atoms the photon is far in the 
wing. The resulting normalized probability distribution 
is 

f{u ^^nH(a,x) (x-u^ + af (13) 
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Since Eq. (| 1 3|) is not analytically integrable, U|| is gen- 
erated fr om this distribut ion by means of the rejection 
method (jPress et al.l 119921 ): a random value of a com- 
parison function that is integrable and lies everywhere 
above f(u\\) is found, and accepted if a second random 
number 6 1Z € [0, 1] (a "univariate" ) is less that the ratio 
of the two functions. Due to the peculiar shape of /(mm ) 
(see Fig. O, following IZtieng fc Miralda-Escudel (|2002f l. 
two comparison functions are used. For the wide range 
of temperatures and frequencies involved we find that a 
satisfactory average acceptance-to-rejection ratio of or- 
der unity is achieved for 

( for < x < 0.2 

u = < x - O.Ola 1 /^ 12 * for 0.2 < x < x cw (a) (14) 
I 4.5 for x > x cw (a) 

as the value Uq separating the two comparison functions. 
Here x cw defines the boundary between the core and the 
wings of the Voigt profile, i.e. where e~ x / \/tt — a/irx 2 . 
The solution to this equation can be approximated as 

x cw (a) = 1.59 - 0.60 log a -0.03 log 2 a. (15) 

When x ~ 0, the photon barely diffuses spatially. Only 
when it has diffused sufficiently far in frequency space 
will it be able to make a large journey in real space. 
The photon may scatter thousands or even hundreds of 
thousands of times before entering the wing of the line 
profile. These scatterings are insignificant in the sense 
that they do not contribute to any important displace- 
ment in neither space nor frequency. Hence, we may 
as well skip them altogether and go directly to the first 
scattering that pushes the photon into the wing. This 
highly efficient accelera tion of the code is achieved fol- 
lowing HETiTai] (|2002f) : if | a; | is less than some critical 
value cc cr it, wj_i,2 is drawn from a truncated Gaussian so 
as to favor fast moving atoms and artificially push the 
photon back in the wing. The resulting rando m velocity 
generator can be written (jDiikstra et al.ll2006t ) as 

1/2 

U -L,l = (^crit - In^l) C0S27T^ 2 ^ 

wi,2 = (x 2 rit - lnT^i) 1 ^ 2 sin27r7?.2, 

where TZ\ and 7?. 2 are two univariates. 

However, the value x cr it is not simply equal to x cw , 
since for a nondense medium, a core scattering can in fact 
be associated with a considerable spatial journey, while 
for clouds of extremely high density even scatterings in 
the inner part of the wing may be neglected. Moreover, 
the exact value of x cr it is actually quite important; this 
acceleration scheme can decrease the computational exe- 
cution time by several orders of magnitude but too high 
values will push the photons unnaturally far out in the 
wings, leading to incorrect results. From the Neufeld 
solution we know that the important parameter is the 
product otq. Correspondingly, we expect x cr ;t to be a 
function of the value of otq in the current cell. Indeed, 
it is found that the value 

_ J for ar < 1 , , 

Xcrit ~ \ 0.02e« lnXQTo for ar > 1, 1 '> 

6 Random numbers i n the interval [0, 1 ] are generated by means 
of the subroutine rani JPress et al.lll992T) . 



where = (0.6,1.2) or (1.4,0.6) for ar < 60 or 

ar > 60, respectively, can be used without affecting the 
emergent spectrum in both various tests (Q and realistic 
situations (Sj6|). Of course, if the photon is already in the 
wing, the proper Gaussian velocity distribution is used, 
i.e. x crit = 0. 

The final frequency Xf of the scattered photon (in the 
reference frame of the fluid element) depends on direc- 
tion in which the photon is scattered, given by the phase 
(probability) function 

W{6) cx 1 + — cos 2 9, (18) 

q 

where is the angle between and the outgoing di- 
rection rif, and R/Q is the degree of polarization for 
90° scattering. For reasons of symmetry, the scatter- 
ing must always be isotropic in the azimuthal direction 
and hence independent of <j). For scattering in the line 
center, i.e. for x < x cw , transitions to the 2Pi/ 2 state 
results in isotropic scattering such that R/Q = 0, while 
the 2P 3 / 2 transition causes some polarization, resulting 
in R/Q = 3/7 (jHamiltonl [l 940( 1. Since the spin multi- 
plicity is 2 J + 1, with J the angular momentum of the 
state, the probability of being excited to the 2P 3 / 2 state 
is twice as large as being excited to the 2P 1 ^ 2 state 7 . For 
scatterings in the wing, polarization fo r it/2 scattering is 
maximal, i.e. R/Q = 1 (|Stenfldll980h . 

The transition between the wing and the core is taken 
to occur at x cw . Obviously, the phase function does not 
change abruptly at this point, but rather varies contin- 
uously from one to another in some fashion. However, 
the difference in the final outcome is very small, even if 
a single phase function is used for all scatterings. 

In the observers frame, the final frequency is then 

x f = Xi - u\\ + hf ■ u + g(l - h t ■ hf), (19) 

where the factor g = hp\Vo/mncvth ((Field! [l959l l. with 
hp\ the Planck constant, accounts for the recoil effect. 

After the scattering, the photon is assigned a new, ran- 
dom value of r and continues its journey. 

3.4. Observations 

Following the above scheme, the photon is trailed as 
it scatters trough real and frequency space, until even- 
tually it escapes the computational box. Subsequently, 
this procedure is repeated for a number (n p h in total) 
of photons sufficiently large that the result converges. 
Here, "convergence" is defined as the change in the de- 
sired result after n p h photons as compared to the re- 
sult after n p h photons being on the ~ 1% level. If we 
merely concern ourselves with the characteristic angu- 
larly averaged spectrum of escaping Lya photons, n v h 
needs not be very large, of the order 10 3 . However, 
since in general the morphology of a galaxy may very 
well cause an anisotropic luminosity, it may be more in- 
teresting to see how the system would appear when ob- 
served from a given angle, at a distance given by the 
redshift of the galaxy. Because the ratio of photons es- 
ca ping in a particular direction i s effectively zero, follow- 
ing l^suFSadi^&^^orrii l|l984f ) we calculate instead for 

7 For the environments produced here, transitions to the 25 
state and subsequent destruction of the photon through two-photon 
processes can be neglected. 
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each scattering and for each photon the probability of es- 
caping the medium in the direction of the observer, or, 
in fact, six different observers situated in the positive 
and negative directions of the three principal 
W(9)e~ Tasa , where 9 is now given by the angle between 
fii and the direction of the observer and T esc is the optical 
depth of the gas lying between the scattering event and 
the edge of the computational box (integrated through 
the intervening cells and of course taking into account 
the different bulk velocities of the cells). 

This probability is added as a weight to a three- 
dimensional array ( "CCD" ) of two spatial and one spec- 
tral dimension. Each pixel suspends a solid angle O p i X 
of the computational box. The total surface brightness 
SBpi x of the area covered by the pixel, measured in en- 
ergy per unit time, per unit area at the location of the 
observer, per unit solid angle suspended by the pixel is 
then 

Ltnt 1 



(20) 



ph., scat. 

where is the luminosity distance given by the redshift , 
and the sum is over all photons and all scatterings. Note 
that Eq. ([20]) does not contain a factor l/4ir, due to the 
fact that the phase functions are normalized to unity. 

Eq. (f20|) is the SB that an observer would measure at 
a distance from the galaxy. Hence, this is the inter- 
esting quantity for comparing with actual observations. 
Theorists tend to be more concerned with the intrinsic 
SB, i.e. the flux measured by a hypothetical observer at 
the location of the source. In this case Eq. (f20|) must 
be multiplied by a factor (1 + z) 4 , and the SB is then 
measured in energy per time per area. 

When a sufficient number of photons has been propa- 
gated, the 3D array can be collapsed along the frequential 
direction to give a "bolometric" Lya SB map, along the 
two spatial directions to give the integrated spectrum, or 
along all directions to give the total flux received from 
the source. Since in fact we obtain a full spectrum for 
each pixel, it is also possible to perform 2D (long slit) 
spectroscopy, giving frequency as a function of position 
of a selected part of the image. 

For the spectra to converge, n p h should be of order 10 4 . 
The luminous regions of the SB maps need n p h ~ 10 5 
to converge, while the outer regions need several 10 6 . 
However, when smoothing the maps so as to simulate the 
effect of the atmosphere, or in any case a finite angular 
resolution, and averaging the SB maps in the azimuthal 
direction to produce SB profiles, less photons are needed, 
of order 10 5 . In the simulations described in $6l n p h ~ 
10 6 -10 7 was used. 

4. TESTING THE CODE 

4.1. Individual Scatterings 

The various probability distribution generators were 
tested against their analytical solutions (in the case such 
solutions exist; otherwise against numerical integration). 
We show here only the result for the parallel velocities 
"II (Fig.©. 

To test the individual scatterings, Fig. [3] shows the re- 
lation between the frequency of the incident and of the 
scattered photon, compar ed with the exac t redistribution 
function as formulated bv lHummerl (|1962f) . Furthermore, 
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Fig. 2. — Probability distribution P(ii||) of parallel velocities 
it 1 1 of the scattering atom for four values x of the frequency of 
the incoming photon, as generated from Eq. JT3}. For photons 
in the line center (x = 0, blue), P(u^) resembles the natural line 
broadening function. For successively larger, but relatively low 
frequencies (x = 2, green, and x = 3, yellow), the photon still has 
a fair chance of being scattered by an atom to which it appears close 
to resonance. For larger frequencies (x = 4.5, red), however, atoms 
with sufficiently large velocities are so rare that the distribution 
instead resembles a regular Gaussian, slightly shifted toward un = 
x. The method for generating Ui is quite good at resolving the 
resonance peak. This is particularly visible in the logarithmic plot 
(bottom panel). 

the rms and mean shift, and the average number of scat- 
tering before returning to the core for wing photons are 
shown. For x — > oo, t he values are seen to converge to 
the results derived by lOsterbrockl (|1962f> . and given by 
Eqs. © and flU). 

4.2. Uniform Slab of Gas 

The most basic confirmation of the reliability of the 
code is a test of the Neufeld solution. Hence, a simu- 
lation of a slab (i.e. with the x- and y-dimension set to 
infinity) is run in which the bulk velocity of the elements 
is set to zero, while the temperature and hydrogen den- 
sity are constant in such a way as to give the desired 
line center optical depth To from the center of the slab to 
the surface. Base cells arc refined in arbitrary locations, 
to an arbitrary level of refinement. Also, W{9) = con- 
stant is used and the recoil term in Eq. (|19p is omitted 
to match the assumptions made by Ncufcld. The result 
for different values of To is shown in Fig. [H while the 
result of varying the initial frequency is shown in Fig. [5l 

For the lowest optical depth {tq — 10 5 , corresponding 
to cltq = 47 at T = 10 4 K), the fit is not very accurate. 
However, this is not an artifact caused by, e.g., a too low 
number of photons in the simulation, but merely reflects 
the fact that the Neufeld solution is no longer valid when 
the optical depth becomes too low (at low optical depths, 
the transfer of photons is no longer dominated by wing 
scatterings, where the line profile can be approximated 
by a power law). 

Figure [S] shows the average number of scatterings 
(N sca t)- Of course, in this case a non-accelerated version 
of the code (i.e. x cr n = 0) was used, since we are inter- 
ested in the true number of scatterings. To get a feeling 
for the physical significance of the optical depths, the 
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Fig. 3. — Tests of the relation between the frequency Xi of the incoming photon, and the frequency xj of the outgoing photon (top 
left). For photons close to the line center, frequencies are distributed more or less uniformly over the line profile. For larger x, fr equencies 
close to the incoming frequency are preferred, but also frequencies of opposite sign. The distribution follows that predicted by IHummerl 
(1962). For even larger frequencies, photons are less likely to be scattered by atoms to which they are at resonance, and the outgoing 
frequency is then only a few Doppler widths away from the ingoing. For sufficiently large x, the rms shift (Ax)rms — > 1 (top right), the 



mean shift (Ax ) — > — l/\x\ (bottom left), and the average number of scattering needed to return to the Core -/Vscat,ret. 
as predicted by [Ostcrbrock (196^). 



(bottom right), 




Fig. 4. — Emergent spectrum of photons injected in the line 
center in an isothermal and homogeneous slab of gas, for different 
values of line center optical depth tq from the center of the slab to 
the surface, compared with the corresponding Neufeld solutions. 
For increasingly optically thick media, the photons must diffuse 
in frequency further and further from the line center in order to 
escape the medium. For all simulations, T = 10 4 K (corresponding 
to a = 0.00047) and n p h = 10 was used. The analytical solution 
becomes increasingly more accurate as tq — » oo. 




Fig. 5. — Emergent spectrum of 10 5 photons injected with 
different initial frequencies Xi n j in a slab of line center optical 
depth to = 10 7 and temperature T = 10 4 K (corresponding to 
OTq = 4700), compared with the corresponding Neufeld solutions. 



region of tq is divided into the domains of the so-called 
Lyman-limit systems (LLSs) and damped Lyct systems 
(DLAs), characterized by limiting neutral hydrogen col- 
umn densities of N m = 10 172 cm" 2 and 7V H i = 10 20 3 
cm~ 2 , respectively. 



4.3. Gas Bulk Motion 

Exc ept for slightly different factors, iDiikstra et al.l 
(|2006f) found that the emergent spectrum, its maximum, 
and the average number of scatterings of an isothermal, 
homogeneous sphere with no bulk velocity of the gas re- 
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Fig. 6. — Average number of scatterings {N S c&t) (triangles) for 
different line center optical depths to, compared with the analyt- 
ical solution (red) given by Eq. Q. The dashed lines indicate 
the regions of optical depths for Lyman-limit systems (LLSs) and 
damped Lyman alpha systems (DLAs). A temperature of T = 10 
K was used for all simulations, and the number of photons per sim- 
ulation varied from 10 3 to ~ 10 5 for the highest and lowest optical 
depths, respectively. 

semble those of the slab. Although not shown here, such 
simulations were also carried out, and were run in two 
different ways; with a version of the code that has con- 
centric shells instead of cells, thus exploiting fully the 
spherical symmetry, and a normal, cell-based version, the 
output of which converges to the former for sufficiently 
high resolution. 

To test if the implementation of the bulk velocity 
scheme produces reliable results, we inspect the emergent 
spectrum of a sphere subjected to isotropic, homologous 
expansion or collapse. Thus, the velocity Vb u ik(r) of a 
fluid element at a distance r from the center is set to 



Vbuik(r) = Hr : 



(21) 



where the Hubble-like parameter Ti is fixed such that 
the velocity increases linearly from in the center to a 
maximal absolute velocity u max at the edge of the sphere 
(r = R): 



H 



R ' 



(22) 



with u ma x positive (negative) for an expanding (collaps- 
ing) sphere. 

For T K, no analytical solution for the spectrum 
exists. Qualitatively, we expect an expansion to cause 
a suppression of the blue wing and an enhancement of 
the red wing of the spectrum. The reason for this is that 
photons blueward of the line center that would otherwise 
escape the medium, are shifted into resonance in the ref- 
erence frame of atom lying closer to the edge, while red 
photons escape even more easily. Conversely, a collaps- 
ing sphere will exhibit an enhanced blue wing and a sup- 
pressed red wing. This is indeed seen in Fig. [7J Another 
way to interpret this effect is that photons escaping an 
expanding cloud are, on the average, doing work on the 
gas, thus losing energy, and vice versa for a collapsing 
cloud. 

In Fig. results for a sphere of gas expanding at dif- 
ferent velocities are shown. For increasing ti max , the 
position of the red peak is progressively enhanced and 
displaced redward of the line center. However, above a 



certain threshold value the velocity gradient becomes so 
large as to render the medium optically thin and allow 
less redshifted photons to escape, making the peak move 
back toward the line center again. 

The results matches closely those found by previous au- 
thors (IZheng fc Miralda-Escudel l2002l: lTasitsiomiir2006al : 
IVerhamme et al.ll2006f )~ 

5. SEMIANALYTICAL ACCELERATION SCHEME 

Most of the computing time is spent in the very dense 
cells. Since each cell is in fact a "uniform" cube, i.e. a 
cube of homogeneous and isothermal gas, if an analyti- 
cal Ncufeld-cquivalent solution for the distribution of fre- 
quency exists, it would be possible to skip a great number 
of scatterings and thus speed up the code further. 

The slab solution is an alternate series which can be 
written in closed form. Unfortunately, this is not feasible 
for the cube soluti on, but under certain approximations, 
iTasitsioinil ()2006b[ ) found that it is still possible to write 
it as an alternate series. The problem is that, whereas 
for the slab the terms quickly die off, the same is not 
true for the cube. In fact she found that to achieve an 
accuracy better than 3%, one must exceed 30 terms. 

Hence, it seems more convenient to seek a "Neufcld- 
based" approximation. Since for the cube, the radiation 
can escape from six faces rather than just two, we may 
expect the emergent radiation to be described by a func- 
tion similar to the slab solution, but using a lower value 
of aTQ . 

5.1. Emergent Spectrum 

Toward these ends, a series of simulations is run in 
which photons arc emitted isotropically from the center 
of a cube of constant — but different — temperature 
and density, and zero bulk velocity. The distance from 
the center to each face is zq. We will investigate opti- 
cal depths To = 10 5 , 10 6 , 10 7 , and 10 8 (measured along 
the shortest path from center to face). In all simula- 
tions, riph = 10 5 , and different temperatures are tested. 
We then fit a Neufeld profile to the emergent spectrum, 
using r\a,TQ as the independent variable, where r\ is the pa- 
rameter to be determined. A priori, we have no reason 
to believe that the same value of 77, if any, should be able 
to describe all optical depths. However, it is found that, 
save for the lowest optical depth (tq ~ 10 5 ), excellent fits 
are obtained using 

r] = 0.71. (23) 

This is seen in Fig. O 



5.2. Directionality of the Emergent Photons 

In realistic, cosmological simulations, the direction 
with which the photons exit the cell is also important. 
Since in the limit to -t 00, any finite size step not 
perpendicular to the surface will just shift to position 
of the photons in the parallel direction, for extremely 
optically thick slabs, the photons should have a ten- 
dency to exit perpendicula r to the surface. In this case, 
iPhillips fc Meszarosl (jl986t ) found that the directionality 
of the emergent radiation approaches that of Thomson 
scattered radiation from electrons, with intensity 



(24) 
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Fig. 7. — Emergent spectrum from an isothermal (T = 10 4 K) and homogeneous sphere of gas undergoing isotropic expansion (red) or 
contraction (blue) in such a way that the velocity at the edge of the sphere is ri ma x ±200 km s _1 . Left panel shows the result for a column 
density -/Vjj : from the center to the edge of 2 X 10 18 cm -2 , corresponding to tq = 1.2 X 10 5 and characteristic of a typical LLS. Right panel 
shows the result for ./Vjjj = 2 X 10 20 cm -2 (to = 1.2 X 10 7 ), characteristic of a typical DLA. Als o shown is the r esult from a simulation 
with iibuik = (black dashed), and the analytical solution for the static sphere (green) as given by [Dijkstra ct al. (2006). For the LLS, to 
is clearly too small to give an accurate fit. 

where \i = cos 9, with the angle between the outgoing 
direction n/ of the photon and the normal to the surface. 

Since the number of photons emerging at fj, is oc 
I{fi)fjL d[i, the probability P( < fi) of exiting the slab with 
fi < fi' is (|Tasitsiomill2006bh 
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Fig. 8. — Emergent spectrum from an isothermal (T = 10 4 K) 
and homogeneous sphere of hydrogen column density ./Vjjj = 2 X 
10 20 cm -2 (a DLA) undergoing isotropic expansion with different 
maximal velocities » max at the edge of the sphere. For increasing 
Umax i the peak of the profile is pushed further away from the line 
center. However, if i; ma x becomes too large, the medium becomes 
optically thin and the peak moves back toward the center again. 




Fig. 9. — Emergent spectra of a uniform cube of different optical 
depths. Neufcld profiles arc fitted to the spectra using ??aro, with 
rj = 0.71 for all tq. 



J M (1 + 2/x)/^ 



^(3 + V). 



(25) 



We confirm that this is also an excellent description for 
a cube (Fig. [TO]). 

The probability distribution is found by differentiating 
Eq. (|25p and recognizing that (i must be positive for the 
photon to escape: 



f(/i + 2/! 2 ) for < n < 1 
otherwise, 



(26) 



Since Eq. (|26|) is valid for all six faces of the cube, the 
azimuthal angle (j> parallel to the face canno t, as in the 
case of a slab, be evenly distributed in [0, 2ir] (jTasitsiomil 
I2006bl) . However, as can be seen from Fig. [TOl the devia- 
tion from uniformity is quite small, and can probably be 
neglected. Furthermore, it seems less pronounced, the 
higher the optical depth. 

5.3. Point of Escape 

The final parameter characterizing the photons escap- 
ing the cube is the point x esc where it crosses the face. 
Figure fTTI shows the azimuthally averaged SB profiles of 
the emergent radiation as a function of distance from the 
center of the face, for different optical depths. It is found 
that the SB profile is fairly well described by a truncated 
Gaussian 

e -(r/z f/2a 2 SB f or o < r < z V2 



SB(r/z ) = < V2^sb 



for r > ZoV2. 



(27) 
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Fig. 10. — Directionality of the photons emerging from a unifo rm cube for different values of to. In the direction perpendicular to the 
face of the cube (left), n j follows the distribution given by Eq. (12(it . while in the azimuthal direction (right) there is a slight deviation 
from isotropy. 

face of the host cell), the normal scheme is used. Other- 
wise, the photon is assigned a new frequency according to 
the effective Neufeld distribution: drawing a univariate 
TZ and setting this equal to the Neufeld-cquivalent cube 
solution 8 integrated from -co to i yields (after some 
algebra) 




Jcube(T ,x)dx 



= _ tan -1 e V w V54(4-*f)MaT )eff 



(28) 



Inverting the above expression, the frequency x of the 
photon then becomes 



Fig. 11. — Probability distribution (solid lines) of the exiting 
point for photons emerging from a uniform cube of side length 2zo, 
as a function of distance r from the center of the face, normalized 
to zo , for different values of to . The distributions have been calcu- 
lated as best fits to the corresponding simulated surface brightness 
profiles (dotted lines), as given by Eq. 1271 . 

The dispersion ctsb of the SB decreases very 
slowly with optical depth, and can be written as 
ctsb = 0.48 - 0.041ogT /10 8 . However, in the con- 
text of a cell-based structure, one might state that 
it is meaningless to discuss differences in position on 
scales smaller than the size of a cell, and it is found that 
final results are not altered by simply setting ctsb = 0.5. 



5.4. Implementation of the Cube Solution 

With the probability distributions of frequency, direc- 
tion and position for the photons escaping the cell, we 
are now able to accelerate the code further: every time 
a photon finds itself in a host cell of aro higher than 
some given threshold, which to be conservative we define 
a s £J7"o > 2 x 10 3 , an effective cell with the photon in 
the center is built, with "radius" zq equal to the distance 
from the photon to the nearest face of the host cell. Since 
the effective cell is always completely circumscribed by 
the host cell, its physical parameters are equal to those 
of its host cell. 

If the value of aro in the effective cell, (aro) e ff, is be- 
low the threshold (i.e. if the photon is too close to the 



x f 




54 . . , nil 
VKO-TojeS m tan — 



1/3 



(29) 



Finally, the direction and the position of the photon is 
determined from Eqs. (j2l)|) and (|27[) (with a probability 
of escaping from a given face equal to 1/6), whereafter it 
continues its journey. 

6. APPLICATION TO COSMOLOGICAL SIMULATIONS 

To demonstrate the potency of the developed code, 
we apply it to a number of different simulated galax- 
ies. Specifically, we will study three young "Lyman-break 
galaxies" (LBGs), for the purpose of the present study 
dubbed K15, K33, and S115, at a redshift of z = 3.6, at 
which time the Universe was ~ 1.8 Gyr old. These three 
galaxies are representative of typical galaxies in the sense 
that the two larger evolve into Milky Way/M31-like disk 
galaxies at z = 0, whereas S115 becomes a somewhat 
smaller disk galaxy. The characteristic circular velocities 
of the galaxies at z = are V c = 245, 180 og 125 km s _1 , 
respectively. 

The cosmological simulation is conducted using an 
N-body/hydrodynamical TreeSPH code. A "pseudo"- 
RT scheme of the ionizing UV radiation is accomplished 
on-the-fly, while a more comprehensive UV RT scheme 
is implemented post-process. This is described below. 



Of course normalized to unity instead of the usual 1/4-7T. 
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6.1. Underlying Cosmological Simulations 

The numerical simulations are first carried out at low 
resolution, but in a large volume of space. Subsequently, 
interesting galaxy forming regions are rcsimulated at 
high resolution. Typically, resimulations are performed 
at 64 x higher mass resolution. In order to check the sig- 
nificance of the resolution, also lower resolution (8x) and 
ultrahigh resolution (512 x) simulations are executed. 

The simulations are started at an initial rcdshift Zi = 
39, at which time there is only dark matter (DM) and 
gas SPH particles. The latter eventually evolves partly 
into star particles, while, in turn, star particles can be- 
come gas pa rticles again. The star forma tion criteria is 
describ ed in iSommer-Larsen et al.l (|2003f ). For K15 and 
K33, a iKroupal (|l998f ) initial mass function (IMF ) has 
been assumed, while for S115, a ISalpeterl (|1955l ) IMF 
was used. A standard ACDM is assumed, i.e. fl m = 0.3, 
Oa = 0.7, and the rms linear density fluctuation on scales 
of 8/i _1 Mpc is <j$ = 0.9. In addition to hydrogen and 
helium, the code also follows the chemical evolution of 
C, N, O Mg, S i, S, Ca, and Fe, using the method of 
iLia et all (f2002h . 

Additional information about the simulations can be 
found in Tab. [TJ while Tab. [2] summarizes the physical 
properties of the resulting galaxies, demonstrating that 
they are typical of galaxies at z = 3-5. Having an R- 
band magnitude of R ~ 25, K15 and K33 could be de- 
tected as LBGs in a survey like the one of ISteidel et alJ 
(120031 R nm = 25.5). With R ~ 28, S115 is too faint 
to be detected as an LBG in cur rent surveys, although 
the li mit is being approached (e.g. ISawicki fe Thompson! 
120051 reaching R Vlm ~ 27). 

The Lya emissio n is produced by three d i fferen t 
processes (see also lLaursen fe Sommer-Larsenl I2007D . 
viz. from recombinations in photoionized regions around 
massive stars (responsible for ~ 90% of the total Lya 
luminosity), gravitational cooling (~ 10%), and a meta- 
galactic UV background (UVB) photoionizing the exter- 
nal parts of the galaxy (~ 1%). In the first case, th e lu- 
minosity is determined foll owing iFardal et alJ (l2001f ). us- 
ing the code Starburst99 (jLeitherer et al.lll999f ) to yield 
the Lyman continuum (LyC), and assuming a mean LyC 
photon energy of 1.4 Rydberg and that 0.68 Lya pho- 
tons are emitted per photonionization. The gravitational 
cooling is accounted for by keeping track of temperature 
and ionization state of the gas. The UVB field is assumed 
to be that given bv lHaardt fe Mad"aul (|1996h . where the 
gas is treated as optically thin to the UV radiation until 
the mean free path of a UV photon at the Lyman limit 
becomes less than 1 kpc, at which point the gas is treated 
as optically thick and the UV field is "switched off" . 

For a more thoro ugh description of th e code , the 
re ader is referred t o ISommer- Larsen et afl (|2003[) . and 
to ISommer-Larser] (|2006l) for recent updates. 

6.2. Ionizing UV Radiative Transfer 

To model the propagation of ionization f ronts realisti- 
cally, iRazoumov fe Sommer-Larsenl (|2006l 120071 ) imple- 
mented the following RT scheme as a post-process to the 
cosmological simulation: the position of the SPH parti- 
cles and their associated physical properties are interpo- 
lated from the 50 nearest neighboring particles onto a 
grid of base resolution 128 3 cells. Dense cells are sub- 



divided in eight cells, which are further refined until no 
cell contains more than ten particles. 

Around each stellar source, a system of 12 x 4™" 1 
radial rays is built that split either as one moves far- 
ther away from the source or as a refined cell is en- 
tered, and n = 1,2,... is the angular resolution level. 
Once a radial ray is refined angularly, it stays refined 
at larger distances from the source, even when leaving 
the high-resolution region. In each cell, the photore- 
action number and energy rates due to photons trav- 
eling along ray segments passing through that cell are 
accumulated. These rates are then used to update tem- 
perature and the ionization state of hydrogen and he- 
lium, which in turn are used to calculate the LyC opac- 
ities used in the RT. In addition to stellar photons, we 
also account for ionization and heating by LyC photons 
originating out side the computational vol ume with the 
FTTE scheme H azoumov fe Cardalll[2005h assuming the 
Haardt-Madau UVB. 

Although the UV RT is not coupled to hydrodynam- 
ics, hydrodynamical (shock) heating needs to be taken 
into account during the RT. All cells with temperatures 
above T cr = 3 x 10 4 K in the SPH output are considered 
to be shock-heated, and during the RT only their ion- 
ization state, not their temperature, is updated. For all 
cells with temperatures below T cr a hydro-heating term 
is computed, which is defined as the amount of heating 
needed to keep the temperature of that cell constant if its 
ionization state stayed at the original level. This hydro- 
heating term is used to update both temperature and 
ionization of all T < T CT cells. For these cells, a tempera- 
ture ceiling of T cr is used to avoid unphysical overheating, 
as during the RT calculation heated gas is not allowed 
to expand. 

6.3. Lya Radiative Transfer 

For the purpose of the Lya RT, the same AMR grid is 
used as for the UV RT. A (250 kpc) 3 box is cut out from 
each z = 3.6 snapshot and is refined to £ max = 7 or 8 
levels, where C = corresponds to the unrefined (base) 
grid. Thus, the cell size of the smallest cells is merely 
~ 10 pc, more than four orders of magnitude smaller 
than the computational box itself and comparable to the 
size of molecular clouds. 



cance of the Improved UV RT 



6.3.1. 

Initially, we will focus on K15, consisting of two rather 
compact "disks" embedded in a more extended, ~ 10- 
15 kpc thick, sheet-like structure composed of nonstar- 
forming H I gas, taken to constitute the xy-plane. Figure 
[T2l shows the general effect of the scattering: the SB is 
increased in the outskirts of the galaxy, at the expense 
of a decrease in the center, where most of the photons 
are produced. The effect on the spectrum is also seen: 
while the only broadening of the input spectrum visible is 
due to the bulk motion of the gas elements emitting the 
photons (the natural broadening being much smaller), 
the scattered spectrum is severely broadened, diminished 
by more than an order of magnitude, and split up into 
two peaks due to the high opacity of the gas for photons 
in the line center. 

Figure [LJ] shows the emergent spectrum, as observed 
when viewing the sheet edge-on and face-on, respectively. 
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TABLE 1 

Characteristic quantities of the simulations 



Galaxy 


K15 




K33 




S115 




Resolution 


8X 


64 x 


8X 


64 x 


64 x 


512x 


ATp.tot 


3.0 x 10 5 


2.2 x 10 6 


1.5 x 10 5 


1.2 x 10 6 


2.1 x 10 5 


1.3 x 10 6 


NSPH 


1.4 x 10 5 


1.0 x 10 6 


7.1 x 10 4 


5.5 x 10 5 


1.0 x 10 5 


6.4 x 10 5 


mSPH,mstar 


7.4 x 10 5 


9.3 x 10 4 


7.4 x 10 5 


9.3 x 10 4 


9.3 x 10 4 


1.1 x 10 4 


"DM 


4.2 x 10 6 


5.2 x 10 5 


4.2 x 10 6 


5.2 x 10 5 


5.2 x 10 5 


6.6 x 10 4 


e SPH,estar 


382 


191 


382 


191 


191 


96 


£ DM 


680 


340 


680 


340 


340 


170 


^min 


20 


10 


20 


10 


10 


5 



Note. — Total number JV p ,tot of particles, number TVgPH of SPH particles only, masses 
m, gravity softening lengths e, and minimum smoothing lengths Z m i n of dark matter (DM), 
gas (SPH), and star particles used in the simulation of the three galaxies K15, K33, and 
S115, for different resolutions. Masses are measured in /i _1 Mq, distances in ft _1 pc. 





Fig. 12. — Results for the galaxy S115 from the Lya scattering code MoCaLaTA. Top left plot shows the Lyo surface brightness (SB) 
map of the emitted radiation in the positive ^-direction, i.e. how the galaxy would look if the photons did not scatter. Almost all the 
photons are emitted in the center, and almost none are emitted further away than 20 kpc from the center. Taking into account scattering 
(top right), the emission is clearly much more extended, while the maximum SB is decreased. This is also seen in the SB profile (bottom 
left), i.e. the azimuthally averaged SB map. Both the true (dotted curves) and the profiles of the SB smoothed with a seeing of 0''8 (solid 
curves) are shown, for both the emitted (blue) and the scattered (red) radiation. The photons scatter not only in real, but also in frequency 
space (bottom right); while the emitted spectrum (blue) is close to a delta function, the escaping spectrum (red) is broadened by many 
Angtroms. Moreover, due to the fact that the hydrogen cross-section is so large for photons in the line center, the spectrum is split up into 
two peaks. 
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TABLE 2 

Physical properties of the simulated galaxies 



Galaxy 


K15 


K33 


S115 


SFR/A/q yr- 1 


16 


13 


0.5 


M*/M 


1.3 x 10 10 


6.5 x 10 9 


2.5 x 10 8 


V c (z = 0)/km s" 1 


245 


180 


125 


r vir /kpc 


47.5 


37.3 


22.0 


Lhya/erg s _1 


3.3 x 10 43 


1.4 x 10 43 


7.0 x 10 41 


^i/.uv/erg s~ x Hz -1 


6.7 x 10 28 


5.5 x 10 28 


3.6 x 10 27 



Note. — Star formation rates (SFRs), stellar masses (M„ ) 
circular velocites (V c ), virial radii (r v j r ), Lya luminosities 
(^Lya)i an d UV luminosities (L^tjv) f° r t ne three simulated 
galaxies K15, K33, and S115. All quoted values correspond 
to a redshift of z = 3.6, except V c which is given for z = 0. 

For each case, four spectra are shown: i) using the tem- 
perature and ionization from the original SPH run, n) 
using the temperature and ionization computed with the 
LyC transfer scheme from Scc. l6.2l applicd to stellar pho- 
tons only, ill) using the temperature and ionization from 
a model in which the UVB is traced with the FTTE 
scheme and there are no stellar photons, and iv) using 
distributions from a model in which LyC RT is computed 
for both stellar and UVB photons. 

As expected, fewer photons escape in the plane of the 
sheet, as the Lya optical depth in this direction is greater 
and more photons scatter and eventually escape in the 
directions perpendicular to the sheet. Note that both 
edge-on and face-on spectra are sensitive to the changes 
in ionization of predominantly low-density regions which 
is computed with LyC radiative transfer. As discussed 
in §2.2| the quantity that determines the shape of the 
spectrum is the product oto oc riHi/T (at a fixed physi- 
cal size) . Generally, a higher temperature will also imply 
a lower density of neutral hydrogen, and vice versa, and 
thus we might expect n^i/T to change rapidly to higher 
or lower values for the improved scheme. However, for 
high density cells the change in the ratio n^/T is mi- 
nor when invoking the improved scheme, in most cases 
of order unity. Only in low density cells is this ratio 
considerably altered, but since > 90% of all scatterings 
take place in high-density cells, the overall effect is small. 
The only notable difference is seen in the inner part of 
the spectrum, which is exactly the part that is created 
by the low-density regions, since photons near the line 
center cannot escape from high-density regions. 

6.3.2. Characteristics of the Emergent Spectrum 

The double peak profile seen in Fig. [13] is characteristic 
of Lya emission lines; the high opacity for photons near 
the line center makes diffusion to the either side necessary 
in order to escape the galaxy. Nonetheless, unlike in 
previous simple models, the intensity in the line center is 
not zero. The photons that contribute to this intensity 
are those produced mainly by gravitational cooling, in 
the outskirts of the systems. 

Figures [14] and [15] display the spectra emerging from 
galaxies K33 and S115 at z = 3.6, in six different di- 
rections. The exact shape varies quite a lot, but all 
spectra appear to exhibit the double peak profile. More- 
over, they are broadened by several 1000 km s _1 . 

Double peaks have been observed on several occasions 



(e.g.. lYee fc De Robertisl 119911: IVenemans et~aLl |2005|) . 
iTapken et all (|2007l ). using a resolution of R ~ 2000, 
found three out of 16 LAEs at redshifts z ~ 3-4 
to have a double-peaked profile, while Yamada et al. 
(priv. comm.), using R ~ 1500 found that 26 of 94 LAEs 
at redshift 3.1 have the double-peaked profile. The dif- 
ference in magnitude of the two peaks can be a signature 
of infalling/outflowing gas, cf. £|4.3t In principle, this dif- 
ference may be used as a probe of the gas dynamics, and 
has indeed been used to infer the presence of galactic 
superwinds. However, since LAEs are often situated in 
highcr-than-average density regions, the removal or di- 
minishing of the blue peak might also be caused by IGM 
resonant scattering combined with cosmic expansion. 

Even if the double peak survives intergalactic trans- 
mission, fairly high resolution is required. With a typ- 
ical separation AA of the peaks of the simulated spec- 
tra from a few to ~ 10 A, the resolution must be 
R = 5600/ A A ~ 500-1500. 

It is interesting that while the Lya profile of the sub- 
Milky Way galaxy K15 also shows a moderate gas in- 
fall, we see some outflow signatures in the spectrum 
of the lower-mass galaxy SI 15 from negative x- and y- 
directions. Fitting a Ncufcld, or Dijkstra, profile to the 
observed spectra can give us an idea of the intrinsic prop- 
erties of the system. Unfortunately, due to the degener- 
acy between column density and temperature, one would 
have to gain knowledge either of the parameters by other 
means to constrain the other (e.g. by inferring column 
density from the spectrum of a coincident background 
quasar, or by assuming a temperature of, say, 10 4 K, 
representative of most of the Lya emitting gas). How- 
ever, dumpiness of the ISM will lower the effective op- 
tical depth, making any inferred value of aro a lower 
bound. This is exactly the reason we need realistic mod- 
els; galaxies are not isothermal, homogeneous slabs. 

6.3.3. Surface Brightness 

Figure [R)] displays the Lya SB maps, integrated over 
the line, of two of the galaxies, K15 and K33. As men- 
tioned already, K15 is embedded in a sheet-like structure, 
lying at the intersection of three filaments of gas. Since 
the bulk of the photons are produced in the central, star- 
forming regions, the total optical depth is larger in the 
direction parallel to the sheet than perpendicular to it, 
and hence we would expect the photons to escape more 
easily in the face-on direction. Similarly, K33 is situated 
in a filament of gas, taken to lie along the z-axis. Here, 
we would expect the photons to escape more easily in the 
x- and y-directions. 

Averaging the SB in the azimuthal direction, the SB 
profiles of the three galaxies, each as viewed from two 
different directions, are shown in Fig. [T7] while Tab. 
summarizes the observed maximum surface brightnesses, 
SB max . In fact, regarding K33, SB maxa: _ turns out 

to be smaller than both SB maXiZ _ and SB maXiZ+ , due to 
the presence in the line of sight of hydrogen clouds with 
little star formation causing a shadowing effect. As seen 
in Tab. [3] the observed SB of a given galaxy varies with 
viewing angle by approximately an order of magnitude. 

This result is intriguing in relation to the classification 
of galaxies. Galaxies are commonly annotated accord- 
ing to the method by which they are selected, and one 
of the mysteries in the context of galaxy formation and 
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Fig. 13. — Emergent spectrum of the galaxy K15, as seen when observing the sheet-like structure in which the galaxy is embedded 
edge-on (left) and face-on (right). Blue lines show the spectrum for the model without the improved UV RT, while green, yellow, and red 
lines show the spectrum when treating the UV RT properly for the stellar sources only, the UV background only, and both, respectively. 
The only real difference is seen in the blue peak of the spectra, which is a bit higher for the improved models. In particular, all models 
indicate a moderate net infall of gas, enhancing the blue-to-red peak ratio. This figure, as well as all following figures, refers to z = 3.6. 
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FlG. 14. — Spectral distribution of the photons escaping the galaxy K33 in six different directions; along the positive (+) and negative 
(— ) x-, y-, and z-direction. The dashed line in the middle of each plot indicates the line center. The lower abscissa gives the redshiftcd 
wavelength of the photons while on the upper abscissa, the wavelength distance from the line center is translated into recessional velocity. 
The resonant scattering of hya is seen to broaden the line by several thousands of km s — 1 . 
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Fig. 15. — Same as Fig. [14] but for the galaxy S115. 
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Fig. 16. — Lyct surface brightness maps (integrated over the line) of the two galaxies K33 (left) and K15 (right) as viewed from the 
negative y-direction and the negative ^-direction, respectively. 
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Fig. 17. — Surface brightness (SB) profiles for the inner 50 kpc of the three galaxies K15 (right), K33 (middle), and S115 (left). For 
each galaxy, the SB is showed as observed from two different angles (blue and red, respectively). Both the "true" SB profiles (dotted) and 
the profiles of the SB convolved with a seeing of 0'.'8 (solid) are shown. Left j/-axis gives the SB that would be seen by an observer at the 
location of the galaxies, while right y-axis gives the SB as observed from Earth. 
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TABLE 3 

Maximum Observed Surface Brightnesses from 
Different Directions 



Galaxy 




K15 




K33 




S115 




log SB max 


,x + 


8.6 x 10- 


-3 


6.2 x 10- 


-3 


1.3 x 10- 


-3 


log SB max 


,X — 


3.7 x 10 


-3 


2.3 X 10 


-3 


1.5 x 10 


-3 


log SB max 


:V + 


9.9 x 10" 


-3 


3.2 X 10 


-2 


7.6 x 10- 


-4 


log SB max 


,y- 


1.2 x 10- 


-2 


2.1 x 10- 


-2 


1.2 x 10- 


-3 


log SB max 




4.9 x 10- 


-2 


3.4 x 10" 


-3 


2.9 X 10 


-4 


log SBmax 




5.3 X 10 


2 


4.7 x 10" 


-3 


4.0 x 10" 


-4 


Ratio 




14.3 




13.6 




5.3 





Note. — Surface brightnesses (SBs, calculated from 
the images convolved with a seeing of / .'8) are measured 
in erg s _1 cm -2 , at the location of the galaxies. The max- 
imally and the minimally observed SBmax's for a given 
galaxy are written in boldface, and the ratios between 
these are given in the lower row. 

evolution is the connection between the different types. 

As already mentioned, sufficiently high column den- 
sities of neutral hydrogen in the line of sight toward a 
bright background source give rise to broad absorption 
lines in their spectra, and may be detected as the so- 
called DLAs. On the other hand, galaxies with high 
enough star formation rates (SFRs) may be detected in 
narrowband searches by an excess of their narrowband 
to continuum flux as LAEs. 

Due to the massive amount of neutral hydrogen, DLAs 
are self-shielded against ionizing radiation and may hence 
be able to cool sufficiently to initiate star formation. 
The assumption that DLAs be progenitors of present- 
day galaxies therefore seems reasonable, and is indeed 
generally accepted. Pursuing this idea, one may thus im- 
age an evolutionar y sequence, for instance DLA — » LAE 
-> LBG (see, e.g.. lGawisenl2006at iGawiser et al.ll2007l ). 
However, high-rcdshift galaxy classification may also be 
a simple consequence of a selection effect, i.e. reflecting 
the means by which they are probed. It seems safe to 
say that the exact relation between the different types 
remains unclear. 

All three galaxies of the present study contain enough 
neutral hydrogen to make them detectabl e as DLAs in 
the s pectra of (hypothetical) quasars (see lEllison et alJ 
|2007| ). More interestingly, the present results show that 
while their high SFRs (K15: 16 M© yr" 1 ; K33: 13 M s 
yr -1 ; SI 15: 1/2 M Q yr -1 ) may make at least the two 
larger galaxies detectable as LAEs when viewed from a 
given direction, it may not be possible to see them in 
Lya from another direction; instead, it may be possible 
to observe them as LBGs. This effect demonstrates how 
galaxies selected by different means may be connected to 
each other, although this obviously has to be quantified 
through detailed modeling, including the effect of dust. 

Overlaps in the properties of LAEs and LBGs have 
also been inferred observationallv TGawiser et al.l (|2006bf ) 
found that more than 80% of a sample of cmission- 
line-selected LAEs have the right UVB-co\ots to be se- 
lected as LBGs. The primary difference between the 
two populations is the selection criteria, as only ~ 10% 
are also brighter than the Rab < 25.5 "spectroscopic" 
LBG magnitud e cut. Also, when correcting for dust, 
iGronwall et all (|2007f) found comparable SFRs for the 



two populations. 

For high-redshift LAEs, SFRs are inferred almost 
exclusively from Lya flux measurements, assuming 
isotropic luminosity. However, as is evident from the 
above discussion, in general the complex morphology of 
a galaxy may very well cause a preferred direction of pho- 
ton escape. For the three galaxies of the present study, 
the angular variation in flux can be as high as a factor 
of 3.4, 6.2, and 3.3, for K15, K33, and S115, respectively 
(here the flux is calculated by integrating the SB maps 
over a region of radius r = 25 kpc, centered at r = 0). Al- 
though not as pronounced as in the case of SB max , this 
introduces a considerable source of uncertainty, which 
may propagate into estimates of SFRs or into calcula- 
tions of the content of dust residing in galaxies. 

6.3.4. Testing Resolution and Interpolation Scheme 

In order to check the impact of the resolution of the 
cosmological simulation on the results of the Lya RT, 
the above RT calculations were also carried out on the 
output of cosmological simulations performed at eight 
times lower resolution. Furthermore, the procedure by 
which the SPH particles were interpolated onto the grid 
was tested by using the 10 nearest neighboring particles 
instead of the usual 50 particles, and running similar 
RT calculations on the output grids. Although in both 
cases the results changed somewhat, there seems to be 
no general trend. Due to the slightly different evolution 
of the lo-res galaxies, the precise configuration of stars 
and gas clouds will not necessarily be the same, and thus 
luminous peaks in the SB maps cannot be expected to 
coincide exactly. However, the maximum SBs appear to 
agree to within a few tens of percents, as do the slopes 
and the overall amplitudes of the SB profiles. The out- 
come of three such simulations can be seen in Figures 
□landCEi 

7. SUMMARY AND DISCUSSION 

A Monte Carlo Lya radiative transfer code has been 
presented, and subsequently tested against various an- 
alytical solutions. The code is capable of propagating 
Lya radiation on an adaptively refined mesh, with an 
arbitrary distribution of Lya source emission, gas tem- 
perature, density, and velocity field, and is thus suitable 
for making predictions about the diffusion of Lya radia- 
tion in simulated galaxies of arbitrarily high resolution. 

A similar non-AMR version of the code has already 
verified the extendedness of Lya surface brightness of 
z ~ 3-4 compared to the con tinuum as being due to the 
reson ant scattering of Lya (jLaursen fc Sommer-Larsenl 
2007). In the present work, we address the question of 
the impact of the viewing angle on the observed prop- 
erties of galaxies. We find that the anisotropic escape 
of photons may cause the maximum observed SB to vary 
quite a lot — in the three studied cases on the average by 
approximately an order of magnitude — while the total 
observed Lya flux varies somewhat less; a factor of 3-6. 
We propose that this effect may sometimes cause confu- 
sion in the classification of galaxies, in the sense that the 
same galaxy could be detected as an LBG from one direc- 
tion and a LAE from another. In addition, this angular 
variation may act as a source of error when inferring star 
formation rates. 
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Fig. 18. — Comparison spectra for the resolution test and the interpolation test. Left panel shows the spectrum escaping in the negative 
y-direction of K33, simulated at high (solid curve) and intermediate (dotted) resolution. Middle panel shows the spectrum escaping in the 
negative x-dircction of S115, simulated at ultrahigh (512x) resolution, but interpolating the physical parameters onto the AMR grid using 
the 50 nearest neighboring particles (solid) and the 10 nearest neighbors (dotted). Right panel shows the spectrum escaping in the negative 
^-direction of K15, simulated at high resolution and interpolating from 50 neighbors (solid), compared with intermediate rcsolution/10 
neighbors (dotted). The differences do not change the results qualitatively. Please note that in the middle (right) panel, the intensity has 
been multiplied (divided) by a factor of 10 in order to use the same scale as for all three galaxies. 
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Fig. 19. — Comparison surface brightness (SB) profiles of the six models from Fig. 1181 at 0"8 (solid and dashed curves) overplotted on 
the true SB profiles (dotted). While performing the simulations at different resolutions may shift some of the luminous regions somewhat 
spatially, the maximum SB and the overall slope remain virtually unaltered. Moreover, modifying the number of neighboring particles used 
for the interpolation scheme seems unimportant. 

galactic medium has not been modeled either. IGM ab- 
sorption may remove the blue peak of the spectrum, since 
the cosmological expansion will eventually shift it into 
resonance of neutral hydrogen in the line of sight, thus 
approximately halving the observed flux. The damp- 
ing wing of the absorption profile may sometimes extend 
into the red wing, attenuating the line even more. On 
the other hand, the effect of IGM may be substantially 
reduced, since the line-of-sight density of Hi absorbers 
decrease with decreasing redshift, and by z = 3.6, the 
Universe is largely ionized. The implementation of dust 
and IGM RT will presented in a forthcoming paper. 



As yet, the developed code includes no effect of dust. 
Considerable amou nts of dust have been inferred to re- 
side in LBGs (e g.. ISawicki fc Yedfl998l: ICalzettil l2001fc 
iTakeuchi fc Ishiill2004HRigopoulou et aljboofl ). Correct- 
ing for dust increases the typical inferred value of the 
SFR by a factor of - 5 dReddv fc Steidell l200l . One 
may naively expect the presence of dust to amplify the 
difference between the SB observed from different angles 
of inclination. Nevertheless, taking Lya transfer into ac- 
count, it is not obvious how observations will be affected 
by dust. In fact, the Lya-to-continuum ratio may even 
be boosted if most of the dust resides in clump together 
with the neutral gas, since Lya will scatte r off the clouds 
while continuum p hotons get absorbed (Neufeld 1991; 
lHansen fc Oh1[20M ). Also, since the presence of dust in- 
duces the formation of hydrogen molecules, it may actu- 
ally lower riHi somewhat, making it easier for the photons 
to escape. 

The transfer of the Lya photons through the inter- 
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